#include <vector>
#include <ext/hash_map>
#include <deque>
#include <numeric>
#include <algorithm>
#include <list>
#include <string>


#define DEBUG_LEVEL 8
#include <debug_func.h>
#include <ultoa.h>
#include <progress_indicator.h>
#include <print_error.h>
//#include <tab_delimited_str.h>
#include <stlhelper.h>
#include <index_string.h>
#include <open_stream_name_with_substitutions.h>
#include <any_element_within_range.h>

//boost
#include <boost/algorithm/string.hpp>
#include <boost/foreach.hpp>


#include "get_arg.h"
#include "genes_sans_orphans.h"
#include "data_input.h"
#include "main.h"

using std::numeric_limits;
using std::make_pair;
using std::pair;
using std::runtime_error;
using std::min;
using std::max;
using std::cout;
using std::endl;
using std::cerr;
using std::cin;
using std::deque;
using std::vector;
using std::map;
using std::list;
using std::string;
using std::sort;
using std::ofstream;
using std::ostream;
using boost::algorithm::replace_all_copy;
using leo_algo::any_element_within_range;


//________________________________________________________________________________________


//  open_species_specific_file


//	    helper function for dump_macro_blocks

//________________________________________________________________________________________
void open_species_specific_file(ofstream& ofs, const string& path, const string& species)
{
	if (!boost::contains(path, "[SPECIES]"))
		throw runtime_error("the parameter path \"" + path + "\" should contain "
								 "\"[SPECIES]\". This will be replaced by the actual "
								 "species names for each file.");

	string file_name = replace_all_copy(path, "[SPECIES]", species);
	ofs.open(file_name.c_str());
	if (!ofs)
		throw runtime_error("Could not open the file \"" + file_name + "\".");
}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	t_micro_block

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
// 	add_ortho_gene_pos
// 
// 		gp1 and gp2 are the two orthologous positions which are adjacent
// 		Both have to be specified because only these out of the entire set 
// 		are orthologous may lead to ambiguities
//  
void t_micro_block::add_ortho_gene_pos(unsigned gp1, unsigned gp2, unsigned ortholog_id)
{
	ortholog_id_end = ortholog_id;
	ortholog_ids.push_back(ortholog_id);
	assert(gp1 != numeric_limits<unsigned>::max());
	assert(gp2 != numeric_limits<unsigned>::max());
	if (gp1 > gp2)
		std::swap(gp1, gp2);
	if (ortho_beg == numeric_limits<unsigned>::max())
	{
		ortho_beg = gp1;
		ortho_end = gp2;
		return;
	}
	ortho_beg = std::min(ortho_beg, gp1);
	ortho_end = std::max(ortho_end, gp2);
}



string t_micro_block::orientation_name() const
{
	switch (orientation)
	{
	case eUNDETERMINED:
		return "undetermined";
	case eFORWARDS:
		return "forwards";
	case eBACKWARDS:
		return "backwards";
	}
	return "";
}


void t_micro_block::dump (ostream& outstrm) const
{
	
	outstrm << "[";
	BOOST_FOREACH(unsigned i, gene_positions) 
	{
		outstrm << i << " ,";
	}
	outstrm << "] " << orientation_name() << endl;
}



void t_micro_block::dump (ostream&			outstrm, 
						  deque<unsigned>&	rel_gene_pos_to_gene_pos) const
{
	outstrm << "rel:{beg:"
			<< ortho_beg
			<< ", end:"
			<< ortho_end << "}" << endl;
	outstrm << "[";
	BOOST_FOREACH(unsigned i, gene_positions) 
	{
		outstrm << i << " ,";
	}
	outstrm << "] " << orientation_name() << endl;
	assert(ortho_beg == numeric_limits<unsigned>::max() ||
		   ortho_beg < rel_gene_pos_to_gene_pos.size());
	assert(ortho_end == numeric_limits<unsigned>::max() ||
		   ortho_end < rel_gene_pos_to_gene_pos.size());
	if (ortho_beg == numeric_limits<unsigned>::max())
		outstrm << "abs:orphan" << endl;
	else
		outstrm << "abs:{beg:"
				<< rel_gene_pos_to_gene_pos[ortho_beg]
				<< ", end:"
				<< rel_gene_pos_to_gene_pos[ortho_end] << "}" << endl;
}



//________________________________________________________________________________________

//	are_orthologs_adjacent

//________________________________________________________________________________________
bool are_orthologs_adjacent(const t_deq_deq_u&	rel_ortholog_sets,
							unsigned			ortholog_id1, 
							unsigned			ortholog_id2)
{
	if (ortholog_id1 == ortholog_id2)
		return true;
	return any_element_within_range(rel_ortholog_sets[ortholog_id1], 
									rel_ortholog_sets[ortholog_id2],
									1U);
}



bool t_micro_block::is_flipped_contiguous(const t_micro_block& other) const
{
	debug_func_cout(5);
	if
			(
				are_orthologs_adjacent(rel_ortholog_sets, 
									   ortholog_id_beg, other.ortholog_id_beg) ||
				are_orthologs_adjacent(rel_ortholog_sets, 
									   ortholog_id_end, other.ortholog_id_beg) ||
				are_orthologs_adjacent(rel_ortholog_sets, 
									   ortholog_id_beg, other.ortholog_id_end) ||
				are_orthologs_adjacent(rel_ortholog_sets, 
									   ortholog_id_end, other.ortholog_id_end)
			)
		return true;

	BOOST_FOREACH(unsigned i, ortholog_ids) 
		BOOST_FOREACH(unsigned j, other.ortholog_ids) 
			if (are_orthologs_adjacent(rel_ortholog_sets,  i,  j))
				return true;
	return false;
	
}













//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	t_macro_block

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
bool t_macro_block::operator && (const t_macro_block& other)
{
	if
			(
				are_orthologs_adjacent(rel_ortholog_sets, 
									   ortholog_id_beg, other.ortholog_id_beg) ||
				are_orthologs_adjacent(rel_ortholog_sets, 
									   ortholog_id_end, other.ortholog_id_beg) ||
				are_orthologs_adjacent(rel_ortholog_sets, 
									   ortholog_id_beg, other.ortholog_id_end) ||
				are_orthologs_adjacent(rel_ortholog_sets, 
									   ortholog_id_end, other.ortholog_id_end)
			)
		return true;
	std::deque<unsigned> my_ortholog_ids;
	BOOST_FOREACH(t_micro_block micro_block, blocks) 
	{
		my_ortholog_ids.push_back(micro_block.ortholog_id_beg);
		my_ortholog_ids.push_back(micro_block.ortholog_id_end);
	}
	std::deque<unsigned> other_ortholog_ids;
	BOOST_FOREACH(t_micro_block micro_block, other.blocks) 
	{
		other_ortholog_ids.push_back(micro_block.ortholog_id_beg);
		other_ortholog_ids.push_back(micro_block.ortholog_id_end);
	}
	BOOST_FOREACH(unsigned i, my_ortholog_ids) 
		BOOST_FOREACH(unsigned j, other_ortholog_ids) 
			if (are_orthologs_adjacent(rel_ortholog_sets,  i,  j))
				return true;
	return false;
}

// merge two macro_blocks together
void t_macro_block::prepend(t_macro_block& other)
{
	assert(ortho_beg != numeric_limits<unsigned>::max());
	assert(ortho_end != numeric_limits<unsigned>::max());
	assert(other.ortho_beg != numeric_limits<unsigned>::max());
	assert(other.ortho_end != numeric_limits<unsigned>::max());
	blocks.splice(blocks.begin(), other.blocks);
	ortho_beg = min(ortho_beg, other.ortho_beg);
	ortho_end = max(ortho_end, other.ortho_end);
	m_size += other.size();
}

// move iter (block) from t_micro_blocks to macro block
void t_macro_block::push_back(t_micro_blocks& micro_blocks, t_micro_blocks::iterator iter)
{
	assert(ortho_end != numeric_limits<unsigned>::max());
	assert(iter->ortho_beg != numeric_limits<unsigned>::max() ||
		   m_size == 0);
	assert(iter->ortho_end != numeric_limits<unsigned>::max() ||
		   m_size == 0);
	
	blocks.splice(blocks.end(), micro_blocks, iter);
	ortho_beg = min(ortho_beg, iter->ortho_beg);
	ortho_end = max(ortho_end, iter->ortho_end);
	m_size += iter->size();
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	t_macro_blocks

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
struct t_macro_blocks : public list<t_macro_block>
{
	// move iter (block) from t_micro_blocks to t_macro_blocks
	void add_to_last_macro_block(t_micro_blocks& blocks, t_micro_blocks::iterator iter)
	{
		debug_func_cout(5);
		back().push_back(blocks, iter);
	}

	// add iter (block) from t_micro_blocks to new t_macro_blocks
	void add_to_new_macro_block(t_micro_blocks& blocks, t_micro_blocks::iterator iter)
	{
		debug_func_cout(5);
		push_back(t_macro_block(blocks, iter));
	}
	bool last_macro_block_flipped_contiguous(const t_micro_block& block)
	{
		return back().is_flipped_contiguous(block);
	}
};

typedef deque<t_macro_blocks> t_macro_blocks_by_chromosome;





//________________________________________________________________________________________

//	construct_macro_blocks

//		merge contiguous but flipped micro blocks into macro blocks

//________________________________________________________________________________________
void construct_macro_blocks(t_macro_blocks_by_chromosome&	chromosomes_macro, 
						    t_micro_blocks_by_chromosome&	chromosomes_micro)
{
	debug_func_cout(5);

	chromosomes_macro.clear();

	BOOST_FOREACH(t_micro_blocks& chromosome_micros, chromosomes_micro) 
	{
		t_micro_blocks::iterator block	= chromosome_micros.begin();
		t_micro_blocks::iterator end	= chromosome_micros.end();

		// add macro block for this chromosome
		chromosomes_macro.push_back(t_macro_blocks());
		t_macro_blocks& chromosome_macro = chromosomes_macro.back();

		// empty chromosome: ignore
		if (block == end)
			continue;


		t_micro_blocks::iterator next = block; ++next;

		// add first micro block to new macro block
		chromosome_macro.add_to_new_macro_block(chromosome_micros, block);
		
		while (next != end)
		{
			block = next++;
			// if contiguous, add to current macro block
			if (chromosome_macro.last_macro_block_flipped_contiguous(*block))
				chromosome_macro.add_to_last_macro_block(chromosome_micros, block);

			// else add to new macro block
			else
				chromosome_macro.add_to_new_macro_block(chromosome_micros, block);
		}
	}
}




























//________________________________________________________________________________________


//  is_FORWARDS_adjacent
// 
//	    helper for is_adjacent

//________________________________________________________________________________________
bool is_directional_adjacent(	const deque<unsigned>&		all_prev_pos, 
								const deque<unsigned>&		all_curr_pos, 
								int							direction,
								unsigned&					rel_pos1, 
								unsigned&					rel_pos2,
								const vector<t_strand>&		rel_gene_pos_to_strand,
								t_strand					my_strand)
{
	debug_func_cout(5);

	//
	// get adjacent genes via round trip to orthologs which preserve directionality
	// 		and strand
	//
	if (rel_gene_pos_to_strand.size())
	{
		BOOST_FOREACH( unsigned prev_pos, all_prev_pos)
			BOOST_FOREACH( unsigned curr_pos, all_curr_pos)
			{
				// genes in different ortholog_sets should not have the same ortholog
				assert(curr_pos != prev_pos);

				if (int(curr_pos) - int(prev_pos) == direction &&
					rel_gene_pos_to_strand[curr_pos] * direction == my_strand)
				{
					rel_pos1 = prev_pos;
					rel_pos2 = curr_pos;
					return true;
				}

				if (curr_pos < prev_pos)
					continue;
				else
					break;
			}
		return false;
	}

	//
	// Ignore strand
	//
	BOOST_FOREACH( unsigned prev_pos, all_prev_pos)
		BOOST_FOREACH( unsigned curr_pos, all_curr_pos)
		{
			if (int(curr_pos) - int(prev_pos) == direction)
			{
				rel_pos1 = prev_pos;
				rel_pos2 = curr_pos;
				return true;
			}

			if (curr_pos < prev_pos)
				continue;
			else
				break;
		}
	return false;
}





//________________________________________________________________________________________

//	is_adjacent

//      see construct_initial_micro_block

//________________________________________________________________________________________
bool is_adjacent(	const deque<unsigned>&		all_prev_pos, 
					const deque<unsigned>&		all_curr_pos, 
					eORIENTATION&				orient, 
					unsigned&					relpos1, 
					unsigned&					relpos2,
					const vector<t_strand>&		rel_gene_pos_to_strand,
					t_strand					my_strand)
{
	debug_func_cout(5);

	// is known direction maintained??
	if (orient == eFORWARDS)
	{
		if(is_directional_adjacent(all_prev_pos, all_curr_pos, 1, relpos1, relpos2, 
								   rel_gene_pos_to_strand, my_strand))
		{
			orient = eFORWARDS;
			return true;
		}
		return false;
	}
	if (orient == eBACKWARDS)
	{
		if(is_directional_adjacent(all_prev_pos, all_curr_pos, -1, relpos1, relpos2, 
								rel_gene_pos_to_strand, my_strand))
		{
			orient = eBACKWARDS;
			return true;
		}
		return false;
	}

	//
	// orientation is unknown
	//
	bool forwards	= false;
	bool backwards	= false;
	BOOST_FOREACH( unsigned prev_pos, all_prev_pos)
	{
		BOOST_FOREACH( unsigned curr_pos, all_curr_pos)
		{
			// genes in different ortholog_sets should not have the same ortholog
			assert(curr_pos != prev_pos);
			if (curr_pos > prev_pos)
			{
				if (curr_pos - prev_pos  == 1						&&
					(!rel_gene_pos_to_strand.size()	||
					  rel_gene_pos_to_strand[curr_pos] == my_strand))
				{	
					debug_cout(8, "unknown to forwards: ");
					relpos1 = curr_pos;
					relpos2 = prev_pos;
					forwards = true;
				}
				break;
			}

			// curr_pos < prev_pos
			if (prev_pos - curr_pos  == 1						&&
				(!rel_gene_pos_to_strand.size()	||
				  rel_gene_pos_to_strand[curr_pos] == my_strand))
			{	
				debug_cout(8, "unknown to backwards: ");
				relpos1 = curr_pos;
				relpos2 = prev_pos;
				backwards = true;
			}
		}
	}
	if (!forwards)
	{
		// neither
		if (!backwards)
			return false;

		// backwards only
		orient = eBACKWARDS;
		debug_cout(8, "is backwards only: ");
		return true;
	}

	// forwards only
	if (!backwards)
	{
		orient = eFORWARDS;
		debug_cout(8, "is forwards only: ");
	}

	// both
	return true;
}

//________________________________________________________________________________________

//	get_direction_from_strand

//________________________________________________________________________________________
eORIENTATION get_direction_from_strand(t_strand					mystrand, 	
									   const deque<unsigned>&	all_rel_ortho_pos, 
									   const vector<t_strand>&	rel_gene_pos_to_strand)
{
	if (!rel_gene_pos_to_strand.size() || !all_rel_ortho_pos.size())
		return eUNDETERMINED;

	bool direction[2];
	direction[0] = 
	direction[1] =  false;
	BOOST_FOREACH( unsigned rel_ortho_pos, all_rel_ortho_pos)
		direction[rel_gene_pos_to_strand[rel_ortho_pos] != mystrand] = true;

	if (!direction[0])
		return eBACKWARDS;
	if (!direction[1])
		return eFORWARDS;
	return eUNDETERMINED;
}
				   
//________________________________________________________________________________________

//	construct_initial_micro_block

//________________________________________________________________________________________
void construct_micro_blocks(const t_deq_deq_u&				gene_pos_by_chromosome,
							const t_deq_deq_u&				rel_ortholog_sets,
							const vector<unsigned >&		gene_pos_to_ortholog_id,
							const vector<char>&				gene_pos_to_orphan,
							t_micro_blocks_by_chromosome&	micro_blocks,
							const vector<t_strand>&			rel_gene_pos_to_strand,
							const vector<t_strand>&			gene_pos_to_strand)
{
	debug_func_cout(5);
	assert(gene_pos_to_orphan.size() == gene_pos_to_ortholog_id.size());
	
	unsigned NO_ORTHO = numeric_limits<unsigned>::max();
	t_micro_blocks chromosome_synteny;
	//
	//	all chromosomes
	//
	BOOST_FOREACH(const deque<unsigned>& chromosome, gene_pos_by_chromosome)
	{
		t_micro_blocks	chromosome_micro_blocks;
		unsigned		prev_ortholog_id = std::numeric_limits<unsigned>::max();

		//
		//	one chromosome
		//
		BOOST_FOREACH( unsigned gene_pos, chromosome)
		{
			
			//	ignore orphans
			assert(gene_pos < gene_pos_to_orphan.size());
			if (gene_pos_to_orphan[gene_pos])
				continue;

			//
			// get strand
			//
			t_strand my_strand = gene_pos_to_strand[gene_pos];
			
			//
			//	Add first gene
			//
			if (!chromosome_micro_blocks.size())
			{
				prev_ortholog_id = gene_pos_to_ortholog_id[gene_pos];
				assert(prev_ortholog_id < rel_ortholog_sets.size());
				eORIENTATION strand = 
							get_direction_from_strand(	my_strand, 	
														rel_ortholog_sets[prev_ortholog_id], 
														rel_gene_pos_to_strand);
				chromosome_micro_blocks.push_back(t_micro_block(prev_ortholog_id, 
																rel_ortholog_sets,
																gene_pos,
																strand));
				continue;
			}

			
			//
			//	get positions
			//
			unsigned curr_ortholog_id = gene_pos_to_ortholog_id[gene_pos];
			t_micro_block& curr_block = chromosome_micro_blocks.back();

			// fellow ortholog: add to current block
			if (curr_ortholog_id == prev_ortholog_id)
			{
				curr_block.gene_positions.push_back(gene_pos);
				continue;
			}

			unsigned relpos1 = NO_ORTHO;
			unsigned relpos2 = NO_ORTHO;

			
			// appears to be adjacent: add to current block
			if (is_adjacent(rel_ortholog_sets[prev_ortholog_id], 
							rel_ortholog_sets[curr_ortholog_id], 
							curr_block.orientation, relpos1, relpos2,
							rel_gene_pos_to_strand,
							my_strand))
			{
				curr_block.gene_positions.push_back(gene_pos);
				curr_block.add_ortho_gene_pos(relpos1, relpos2, curr_ortholog_id);
				prev_ortholog_id = curr_ortholog_id;
				continue;
				
			}

			//
			// non-adjacent add to new block
			//
			else
			{
				eORIENTATION strand = 
							get_direction_from_strand(	my_strand, 	
														rel_ortholog_sets[curr_ortholog_id], 
														rel_gene_pos_to_strand);
				chromosome_micro_blocks.push_back(t_micro_block(curr_ortholog_id, 
																rel_ortholog_sets,
																gene_pos,
																strand));
			}
			prev_ortholog_id = curr_ortholog_id;
			assert(prev_ortholog_id < rel_ortholog_sets.size());
		}
		micro_blocks.push_back(t_micro_blocks());
		micro_blocks.back().swap(chromosome_micro_blocks);
	}

	// clean up orphaned blocks
	BOOST_FOREACH(t_micro_blocks& chrm, micro_blocks) 
		BOOST_FOREACH(t_micro_block& block, chrm) 
		{
			assert(block.ortho_beg == NO_ORTHO ||
				   (block.ortho_end != numeric_limits<unsigned>::max() &&
					block.ortho_beg != numeric_limits<unsigned>::max()));

			if (block.ortho_beg != NO_ORTHO)
				continue;
			unsigned gene_pos		= block.gene_positions[0];
			unsigned ortholog_id	= gene_pos_to_ortholog_id[gene_pos];
			block.ortho_beg			= 
			block.ortho_end			= rel_ortholog_sets[ortholog_id][0];
			sort(block.ortholog_ids.begin(),  block.ortholog_ids.end());
			remove_duplicates(block.ortholog_ids);
		}

}



















//________________________________________________________________________________________

//	eliminate_small_neighbouring_blocks

//       helper function for eliminate_small_blocks

//________________________________________________________________________________________
unsigned eliminate_small_neighbouring_blocks(t_micro_blocks_by_chromosome&	chromosomes, 
											 vector<char>&					gene_pos_to_orphan,
											 unsigned 						size,
											 deque<unsigned>&				eliminated_orphans)
{
	debug_func_cout(5);
	unsigned cnt_deleted = 0;
	BOOST_FOREACH(t_micro_blocks& chromosome, chromosomes) 
	{
		t_micro_blocks::iterator block = chromosome.begin();
		t_micro_blocks::iterator end = chromosome.end();
		while (block != end)
		{
			if (block->size() <= size)
			{
				++cnt_deleted;
				// mark all in block as orphaned
				BOOST_FOREACH(unsigned gene_pos, block->gene_positions) 
				{
					assert(gene_pos < gene_pos_to_orphan.size());
					gene_pos_to_orphan[gene_pos] = true;
					eliminated_orphans.push_back(gene_pos);
				}
			}
			++block;
		}
	}
	return cnt_deleted;
}




//________________________________________________________________________________________

//	eliminate_small_blocks

//		returns number of blocks eliminated

//________________________________________________________________________________________
pair<unsigned, unsigned> 
eliminate_small_blocks(	const t_deq_deq_u			gene_pos_by_chromosome[],
						vector<char>&				gene_pos_to_orphan,
						const t_deq_deq_u&			rel_ortholog_sets,
						const vector<unsigned >&	gene_pos_to_ortholog_id,
						unsigned					block_size,
						deque<unsigned>&			eliminated_orphans1,
						deque<unsigned>&			eliminated_orphans2,
						const vector<t_strand>&		rel_gene_pos_to_strand,
						const vector<t_strand>&		gene_pos_to_strand)
{
	debug_func_cout(5);

	//
	//	construct initial micro blocks
	//
	t_micro_blocks_by_chromosome			micro_blocks1;
	t_micro_blocks_by_chromosome			micro_blocks2;
	construct_micro_blocks(	gene_pos_by_chromosome[0], 
							rel_ortholog_sets,
							gene_pos_to_ortholog_id,
							gene_pos_to_orphan,
							micro_blocks1,
							rel_gene_pos_to_strand,
							gene_pos_to_strand);
	construct_micro_blocks(	gene_pos_by_chromosome[1], 
							rel_ortholog_sets,
							gene_pos_to_ortholog_id,
							gene_pos_to_orphan,
							micro_blocks2,
							rel_gene_pos_to_strand,
							gene_pos_to_strand);

	//
	//	eliminate blocks of 1
	//
	return make_pair(	eliminate_small_neighbouring_blocks(micro_blocks1, 
															gene_pos_to_orphan, 
															block_size,
															eliminated_orphans1),
						eliminate_small_neighbouring_blocks(micro_blocks2, 
															gene_pos_to_orphan, 
															block_size,
															eliminated_orphans2));
}




//________________________________________________________________________________________

//	histogram_macro_blocks

//________________________________________________________________________________________
void histogram_macro_blocks(const string&						file_name_template, 
							const string&						species_name, 
							const t_macro_blocks_by_chromosome&	chromosomes,
							map<unsigned, unsigned>&			sizes)
{
	debug_func_cout(5);
	ofstream os;
	open_species_specific_file(os, file_name_template, species_name);

	BOOST_FOREACH(const t_macro_blocks& chromosome, chromosomes) 
		BOOST_FOREACH( const t_macro_block& block, chromosome)
			sizes[block.size()]++;

	typedef pair<const unsigned,  unsigned> t_pairuu;

	// get total
	double inv_total_count = 0.0;
	BOOST_FOREACH(t_pairuu& i, sizes) 
		inv_total_count += i.first * i.second;
	inv_total_count = 1.0 / inv_total_count;



	// print histogram
	double running_count = 0;
	BOOST_FOREACH(t_pairuu& i, sizes) 
	{
		running_count += i.first * i.second;
		os
			<< i.first << "\t"
			<< (inv_total_count * running_count) << "\t"
			<< i.second << "\n";
	}
}

//________________________________________________________________________________________

//	histogram_micro_blocks

//________________________________________________________________________________________
void histogram_micro_blocks(const string&						file_name_template, 
							const string&						species_name, 
							const t_micro_blocks_by_chromosome&	chromosomes,
							map<unsigned, unsigned>&			sizes)
{
	debug_func_cout(5);
	ofstream os;
	open_species_specific_file(os, file_name_template, species_name);

	
	BOOST_FOREACH(const t_micro_blocks& chromosome, chromosomes) 
		BOOST_FOREACH( const t_micro_block& block, chromosome)
			sizes[block.size()]++;

	typedef pair<const unsigned,  unsigned> t_pairuu;

	// get total
	double inv_total_count = 0.0;
	BOOST_FOREACH(t_pairuu& i, sizes) 
		inv_total_count += i.first * i.second;
	inv_total_count = 1.0 / inv_total_count;


	// print histogram
	double running_count = 0;
	map<unsigned, unsigned>		binned_cummulative;
	BOOST_FOREACH(t_pairuu& i, sizes) 
	{
		binned_cummulative[i.first / 10] += i.first * i.second;
		running_count += i.first * i.second;
		os
			<< i.first << "\t"
			<< (inv_total_count * running_count) << "\t"
			<< i.second << "\n";
	}
	os << "\t\n\t \n";
	for (unsigned i = 0; i < 20; ++i)
		os	<< i << "\t"
			<< binned_cummulative[i] << "\t"
			<< binned_cummulative[i] * inv_total_count <<  "\n";
}




//________________________________________________________________________________________

//	list_micro_blocks

//________________________________________________________________________________________
void list_micro_blocks(	const string&					file_name_template, 
						const string&					species_name, 
						t_micro_blocks_by_chromosome&	chromosomes,
						const vector<unsigned>&			gene_pos_to_gene_id,
						const t_index_string&			gene_names)
{
	debug_func_cout(5);
	string file_name = replace_all_copy(file_name_template, "[SPECIES]", species_name);
	ofstream os(file_name.c_str());

	unsigned micro_block_id = 0;
	map<unsigned, unsigned>	sizes;
	BOOST_FOREACH(t_micro_blocks& chromosome, chromosomes) 
		BOOST_FOREACH( t_micro_block& block, chromosome)
		{
			BOOST_FOREACH( unsigned gene_pos, block.gene_positions)
			{
				assert(gene_pos < gene_pos_to_gene_id.size());
				unsigned gene_id = gene_pos_to_gene_id[gene_pos];
				assert(gene_id < gene_names.size());
				os << micro_block_id << "\t" << gene_names[gene_id] << "\n";
			}
			++micro_block_id;
		}
}
















//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



//		merge Macro blocks



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888







void add_to_cluster_recursively(vector<char>&		included, 
								deque<unsigned>&	cluster, 
								const t_deq_deq_u&	adjacency_links,
								unsigned			start_pos,
								unsigned			orig_start_pos)
{
	debug_func_cout(5);
	assert(start_pos < adjacency_links.size());
	BOOST_FOREACH(unsigned i, adjacency_links[start_pos]) 
	{
		debug_cout(5, "adjacent link to " << i);
		assert(i < included.size());
		if (!included[i] && i >= orig_start_pos)
		{
			included[i] = 1;
			cluster.push_back(i);
			add_to_cluster_recursively(included, cluster, adjacency_links, i, orig_start_pos);
		}
	}
}


void get_contiguous_cluster(unsigned				start_pos, 
							const t_deq_deq_u&		adjacency_links, 
							deque<unsigned>&		contiguous_cluster)
{
	debug_func_cout(5);
	debug_cout(5,  "start_pos=" << start_pos);
	vector<char> non_contiguous(adjacency_links.size(), 0);
	non_contiguous[start_pos] = 1;
	while (1)
	{
		// list of items to include should omit non-contiguous
		vector<char> included(non_contiguous);

		// single linkage cluster and sort
		deque<unsigned>	cluster;
		cluster.push_back(start_pos);
		add_to_cluster_recursively(included, cluster, adjacency_links, start_pos, start_pos);
		assert(cluster.size());
		if (cluster.size() == 1)
		{
			contiguous_cluster.swap(cluster);
			return;
		}

		sort(cluster.begin(),  cluster.end());
		#ifndef NDEBUG
		//for (unsigned i = 0; i < cluster.size(); ++i)
		//	std::cerr << "," << cluster[i];
		//std::cerr << "\n";
		#endif

		// find first item which is not contiguous
		unsigned prev_value = cluster.front();
		deque<unsigned>::iterator iter = cluster.begin();
		while (*iter - prev_value < 2 && iter != cluster.end())
			prev_value = *iter++;

		// all contiguous: return
		if (iter == cluster.end())
		{
			contiguous_cluster.swap(cluster);
			return;
		}

		// ignore non-contiguous elements
		while (iter != cluster.end())
		{
			assert(*iter < non_contiguous.size());
			non_contiguous[*iter++] = 1;
		}
	}

}


void merge_macro_blocks(t_macro_blocks& macro)
{
	debug_func_cout(5);
	// get vector of macro blocks 
	vector<list<t_macro_block>::iterator> macro_blocks;
	list<t_macro_block>::iterator iter = macro.begin();
	while (iter != macro.end())
		macro_blocks.push_back(iter++);
	
	// save list of block x-links
	t_deq_deq_u	adjacency_links(macro_blocks.size(), deque<unsigned>());
	for (unsigned ii = 0; ii < macro_blocks.size(); ++ii)
		for (unsigned jj = ii + 1; jj < macro_blocks.size(); ++jj)
			if(*macro_blocks[ii] && *macro_blocks[jj])
			{
				assert(ii < adjacency_links.size());
				assert(jj < adjacency_links.size());
				adjacency_links[ii].push_back(jj);
				adjacency_links[jj].push_back(ii);
			}

	if (!adjacency_links.size())
		return;
	for (unsigned cur_pos = 0; cur_pos < macro_blocks.size(); ++cur_pos)
	{
		deque<unsigned>			contiguous_cluster;
		get_contiguous_cluster(cur_pos, adjacency_links, contiguous_cluster);

		// merge contiguous blocks
		assert(contiguous_cluster[0] == cur_pos);
		if (contiguous_cluster.size() != 1)
		{
			for (unsigned i = 1; i < contiguous_cluster.size(); ++i)
			{
				assert(contiguous_cluster[i] < macro_blocks.size());
				macro_blocks[cur_pos]->prepend(*macro_blocks[contiguous_cluster[i]]);
				macro.erase(macro_blocks[contiguous_cluster[i]]);
			}
		}
		cur_pos = contiguous_cluster.back();
	}
}



//________________________________________________________________________________________

//	dump_macro_blocks

//________________________________________________________________________________________
void dump_macro_blocks(	t_program_args&					args,
						const string&					species,
						t_macro_blocks_by_chromosome&	chromosomes,
						const t_index_string&			chromosome_names1,
						t_index_string&					chromosome_names2,
						const vector<unsigned>&			gene_pos_to_gene_id,
						const vector<float>&			gene_pos_to_geneposext,
						const vector<unsigned>&			gene_pos_to_chromosome,
						const vector<t_strand>&			gene_pos_to_strand,
						const deque<unsigned>&			rel_gene_pos_to_gene_pos,
						const t_index_string& 			gene_names
						)
{


	ofstream file_chrm;
	open_species_specific_file(file_chrm, args.file_chromosome_macro, species);

	ofstream file_full_macro;
	open_species_specific_file(file_full_macro, args.file_full_macro_synteny, species);

	ofstream file_full_micro;
	open_species_specific_file(file_full_micro, args.file_full_micro_synteny, species);

	ofstream file_macro_list;
	open_species_specific_file(file_macro_list, args.file_macro_synteny, species);
	
	unsigned cnt_genes = 0;
	for (unsigned i = 0; i < chromosomes.size(); ++i)
		BOOST_FOREACH( t_macro_block& macro, chromosomes[i])
			cnt_genes += macro.size();
	args.VERBOSE()  << "\t" << cnt_genes << "\tgenes from " << species << "\n";
	

	//
	//	macros per chromosome histogram
	//
	map<unsigned, unsigned>	cnt_macro_per_chromosome;
	BOOST_FOREACH(t_macro_blocks& chromosome, chromosomes) 
		cnt_macro_per_chromosome[chromosome.size()]++;

	typedef pair<const unsigned,  unsigned> t_hist_item;
	BOOST_FOREACH(t_hist_item& item, cnt_macro_per_chromosome) 
		file_chrm << item.first << "\t" << item.second << "\n";


	unsigned macro_index = 0;
	
	for (unsigned i = 0; i < chromosomes.size(); ++i)
	{
		file_full_macro << "Chromosome " << chromosome_names1[i] << "\n";
		file_full_macro << chromosomes[i].size() << "\tcontiguous blocks\n";
		file_full_micro << "Chromosome " << chromosome_names1[i] << "\n";
		file_full_micro << chromosomes[i].size() << "\tcontiguous blocks\n";
		BOOST_FOREACH( t_macro_block& macro, chromosomes[i])
		{
			unsigned ortho_beg_gene_pos = rel_gene_pos_to_gene_pos[macro.ortho_beg];
			unsigned ortho_end_gene_pos = rel_gene_pos_to_gene_pos[macro.ortho_end];
			assert(ortho_beg_gene_pos < gene_pos_to_chromosome.size());
			unsigned ortho_chromosome_id = gene_pos_to_chromosome[ortho_beg_gene_pos];
			assert(ortho_chromosome_id < chromosome_names2.size());
			file_full_macro	<< "\t"		<< macro.size()
							<< "\t[\t"	<< ortho_beg_gene_pos << "\t" 
										<< ortho_end_gene_pos << "]\t"
							<< "chromosome\t"
							<< chromosome_names2[ortho_chromosome_id] << "\n";
			file_full_micro	<< "\t"		<< macro.size()
							<< "\t[\t"	<< ortho_beg_gene_pos << "\t" 
										<< ortho_end_gene_pos << "]\t"
							<< "chromosome\t"
							<< chromosome_names2[ortho_chromosome_id] << "\n";
			BOOST_FOREACH( t_micro_block& block, macro.blocks)
			{
				file_full_micro << "\t\t" << block.size();
				if (block.orientation == eFORWARDS)
					file_full_micro << "\tForwards\t";
				else if (block.orientation == eBACKWARDS)
					file_full_micro << "\tBackwards\t";
				else
					file_full_micro << "\tIndeterminate\t";
				assert(block.ortho_beg < rel_gene_pos_to_gene_pos.size());
				assert(block.ortho_end < rel_gene_pos_to_gene_pos.size());
				unsigned ortho_beg_gene_pos = rel_gene_pos_to_gene_pos[block.ortho_beg];
				unsigned ortho_end_gene_pos = rel_gene_pos_to_gene_pos[block.ortho_end];
				file_full_micro << "[\t" << ortho_beg_gene_pos << "\t" << ortho_end_gene_pos << "]\n";

				file_full_micro << "\t\t\t";
				BOOST_FOREACH(unsigned i, block.ortholog_ids) 
					file_full_micro << "," << i;
				file_full_micro << "\n";
				file_full_micro << "\t\t\t";
				BOOST_FOREACH(unsigned i, block.gene_positions) 
				{
					// save list of genes for this macro block
					file_macro_list << macro_index << "\t" << gene_names[gene_pos_to_gene_id[i]] << "\n";
						
					file_full_micro << "," << gene_names[gene_pos_to_gene_id[i]];
				}
				file_full_micro << "\n";
				file_full_micro << "\t\t\t";
				BOOST_FOREACH(unsigned i, block.gene_positions) 
					file_full_micro << "," << (gene_pos_to_strand[i] > 0 ? 1 : 0);
				file_full_micro << "\n";
			}

			// increment index of current macro block
			++macro_index;
		}
	}
}















//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


int main (int argc, char *argv[])
{
	// we use cout to print progress so flush per call
//	cout.setf(std::ios_base::unitbuf);

	debug_func_cout(5);
//	std::ios::sync_with_stdio(false);

	t_program_args args;
	try
	{
	if (!process_command_line_options(argc, argv, args))
		return 1;


	t_deq_deq_u			gene_pos_by_chromosome[2];
	t_index_string 		species_names;
	t_index_string 		chromosome_names[2];
	t_index_string 		gene_names;
	vector<char>		gene_pos_to_orphan;
	vector<t_strand>	gene_pos_to_strand;
	vector<unsigned>	gene_pos_to_gene_id;
	vector<float>		gene_pos_to_geneposext;
	vector<unsigned>	gene_pos_to_chromosome;
	vector<unsigned >	gene_pos_to_ortholog_id;
	t_deq_deq_u			ortholog_sets;

	read_gene_loci_data(gene_pos_by_chromosome,
						gene_pos_to_orphan,
						gene_pos_to_strand,
						gene_pos_to_gene_id,
						gene_pos_to_geneposext,
						gene_pos_to_chromosome,
						gene_pos_to_ortholog_id,
						gene_names,
						species_names,
						chromosome_names,
						ortholog_sets,
						args.VERBOSE());


	//
	//	sort orthologs
	//
	BOOST_FOREACH(deque<unsigned>& ortholog_set, ortholog_sets) 
		sort(ortholog_set.begin(),  ortholog_set.end());


	//
	//	set unorphaned gene positions
	//
	args.VERBOSE() << "\tGet gene positions without orphans (those on chromosome 'Un' etc) \n";
	vector<unsigned>	gene_pos_to_rel_gene_pos(gene_pos_to_orphan.size(), 0);
	deque<unsigned>		rel_gene_pos_to_gene_pos;
	t_deq_deq_u			rel_ortholog_sets;
	vector<t_strand>	rel_gene_pos_to_strand;
	get_gene_pos_sans_orphans(gene_pos_by_chromosome,
							  gene_pos_to_rel_gene_pos,
							  rel_gene_pos_to_gene_pos,
							  ortholog_sets,
							  rel_ortholog_sets,
							  gene_pos_to_orphan,
							  gene_pos_to_strand,
							  rel_gene_pos_to_strand);

	//
	//	micros histogram
	//
	map<unsigned, unsigned>		micros_histogram1;
	map<unsigned, unsigned>		micros_histogram2;

	//
	//	eliminate genes of less than block size of 1
	//
	if (args.remove_singletons)
	{
		deque<unsigned> eliminated_orphans1;
		deque<unsigned> eliminated_orphans2;
		vector<t_strand> ignore_strand;
		args.VERBOSE() << "\tEliminate small blocks of " << 1 << "\n";
		eliminate_small_blocks(	gene_pos_by_chromosome, 
								gene_pos_to_orphan, 
								rel_ortholog_sets, 
								gene_pos_to_ortholog_id, 
								1,  
								eliminated_orphans1,
								eliminated_orphans2,
								ignore_strand,
								gene_pos_to_strand);
		micros_histogram1[1] = eliminated_orphans1.size();
		micros_histogram2[1] = eliminated_orphans2.size();
		args.VERBOSE() << "\t" << eliminated_orphans1.size() << "\torphaned singles in " << species_names[0] << "\n";
		args.VERBOSE() << "\t" << eliminated_orphans2.size() << "\torphaned singles in " << species_names[1] << "\n";
	}
	

	//
	//	save orphans
	//
	for (unsigned species_id = 0; species_id < 2; ++species_id)
	{
		unsigned orphaned_micro_block_ids = gene_pos_to_gene_id.size();
		ofstream ofs;
		open_species_specific_file(ofs, args.file_orphans, species_names[species_id]);
		BOOST_FOREACH(const deque<unsigned>& chromosome, gene_pos_by_chromosome[species_id])
		{
			BOOST_FOREACH(unsigned gene_pos, chromosome)
			if (gene_pos_to_orphan[gene_pos])
				ofs	<< orphaned_micro_block_ids++ << "\t"
					<< gene_names[gene_pos_to_gene_id[gene_pos]] << "\n";
		}
			
	}

	//
	//	must call this whenever number of orphans changes
	//
	args.VERBOSE() << "\tGet gene positions without orphans and eliminated genes.\n";
	get_gene_pos_sans_orphans(gene_pos_by_chromosome,
							  gene_pos_to_rel_gene_pos,
							  rel_gene_pos_to_gene_pos,
							  ortholog_sets,
							  rel_ortholog_sets,
							  gene_pos_to_orphan,
							  gene_pos_to_strand,
							  rel_gene_pos_to_strand);



	//
	//	construct initial micro blocks
	//
	args.VERBOSE() << "\tInitial micro blocks\n";
	t_micro_blocks_by_chromosome	micro_blocks1;
	t_micro_blocks_by_chromosome	micro_blocks2;
	construct_micro_blocks(	gene_pos_by_chromosome[0], 
							rel_ortholog_sets,
							gene_pos_to_ortholog_id,
							gene_pos_to_orphan,
							micro_blocks1,
							rel_gene_pos_to_strand,
							gene_pos_to_strand);
	construct_micro_blocks(	gene_pos_by_chromosome[1], 
							rel_ortholog_sets,
							gene_pos_to_ortholog_id,
							gene_pos_to_orphan,
							micro_blocks2,
							rel_gene_pos_to_strand,
							gene_pos_to_strand);
	assert(micro_blocks1.size() == chromosome_names[0].size());
	assert(micro_blocks2.size() == chromosome_names[1].size());

	args.VERBOSE() << "\tDump micro blocks\n";
	args.VERBOSE() << "\t" << species_names[0] << " to " << species_names[1] << "\n";
	histogram_micro_blocks(args.file_micro_synteny_histogram, 
				   species_names.string_at(0), micro_blocks1, micros_histogram1);
	list_micro_blocks(args.file_micro_synteny, species_names[0], 
					  micro_blocks1, gene_pos_to_gene_id, gene_names);
	#ifndef NDEBUG
	// dump micro blocks
	//BOOST_FOREACH(t_micro_blocks& chrm, micro_blocks1) 
	//	BOOST_FOREACH(t_micro_block& block, chrm) 
	//		block.dump(cout, rel_gene_pos_to_gene_pos);
	//BOOST_FOREACH(t_micro_blocks& chrm, micro_blocks2) 
	//	BOOST_FOREACH(t_micro_block& block, chrm) 
	//		block.dump(cout, rel_gene_pos_to_gene_pos);
	#endif


	args.VERBOSE() << "\t" << species_names[1] << " to " << species_names[0] << "\n";
	histogram_micro_blocks(args.file_micro_synteny_histogram, 
				   species_names.string_at(1), micro_blocks2, micros_histogram2);
	list_micro_blocks(args.file_micro_synteny, species_names[1], 
					  micro_blocks2, gene_pos_to_gene_id, gene_names);


	//
	//	construct macro blocks which are joined micro blocks
	//
	t_macro_blocks_by_chromosome			macro_blocks1;
	t_macro_blocks_by_chromosome			macro_blocks2;
	args.VERBOSE() << "\tConstruct macro blocks for "  << species_names[0] << endl;
	construct_macro_blocks(macro_blocks1, micro_blocks1);
	args.VERBOSE() << "\tconstruct macro blocks for "  << species_names[1] << endl;
	construct_macro_blocks(macro_blocks2, micro_blocks2);
	assert(macro_blocks1.size() == chromosome_names[0].size());
	assert(macro_blocks2.size() == chromosome_names[1].size());

//	dump_macro_blocks(macro_blocks1, chromosome_names[0]);
//	dump_macro_blocks(macro_blocks2, chromosome_names[1]);
	args.VERBOSE() << "\tMerge macro blocks for "  << species_names[0] << endl;
	BOOST_FOREACH(t_macro_blocks& b, macro_blocks1) 
		merge_macro_blocks(b);
	args.VERBOSE() << "\tMerge macro blocks for "  << species_names[1] << endl;
	BOOST_FOREACH(t_macro_blocks& b, macro_blocks2) 
		merge_macro_blocks(b);
	assert(macro_blocks1.size() == chromosome_names[0].size());
	assert(macro_blocks2.size() == chromosome_names[1].size());

	args.VERBOSE()	<< "\tDump macro blocks for "
					<< species_names.string_at(0) << " to "
					<< species_names.string_at(1) << "\n";

	{
		map<unsigned, unsigned>			sizes;
		histogram_macro_blocks(args.file_macro_synteny_histogram, 
							   species_names.string_at(0), macro_blocks1, sizes);
	}
	{
		map<unsigned, unsigned>			sizes;
		histogram_macro_blocks(args.file_macro_synteny_histogram, 
							   species_names.string_at(1), macro_blocks2, sizes);
	}
	dump_macro_blocks(args, 
					 species_names.string_at(0),  
					 macro_blocks1, 
					 chromosome_names[0], 
					 chromosome_names[1],
					 gene_pos_to_gene_id,      
					 gene_pos_to_geneposext,   
					 gene_pos_to_chromosome,  
					 gene_pos_to_strand,
					 rel_gene_pos_to_gene_pos, 
					 gene_names);
	args.VERBOSE()	<< "\tDump macro blocks for "
					<< species_names.string_at(1) << " to "
					<< species_names.string_at(0) << "\n";
	dump_macro_blocks(args, 
					 species_names.string_at(1),  
					 macro_blocks2, 
					 chromosome_names[1], 
					 chromosome_names[0],
					 gene_pos_to_gene_id,      
					 gene_pos_to_geneposext,   
					 gene_pos_to_chromosome,   
					 gene_pos_to_strand,
					 rel_gene_pos_to_gene_pos, 
					 gene_names);


	args.VERBOSE() << "\n\tAll finished!!\n\n";
	}
	catch(std::exception& e)
	{
		std_print_error(args.CERR(),  e.what());
		return false;
	}
	return 0;
}


/*
//________________________________________________________________________________________

//	dump_micro_block

//________________________________________________________________________________________
void dump_micro_block(string& file_name, t_micro_blocks_by_chromosome& chromosomes, t_index_string& chromosome_names)
{
	ofstream os(file_name.c_str());
	assert(chromosomes.size() == chromosome_names.size());
	for (unsigned i = 0; i < chromosomes.size(); ++i)
	{
		os << "Chromosome " << chromosome_names.string_at(i) << "\n";
		BOOST_FOREACH( t_micro_block& block, chromosomes[i])
		{
			os << block.size();
			if (block.orientation == eFORWARDS)
				os << "\tForwards\n";
			else if (block.orientation == eBACKWARDS)
				os << "\tBackwards\n";
			else
				os << "\tIndeterminate\n";
	
		}
	}
}

//________________________________________________________________________________________

//	are_blocks_contiguous

//________________________________________________________________________________________
inline bool are_blocks_contiguous(const t_micro_block&	b1, 
								   const t_micro_block&	b2)
{
	debug_func_cout(5);
	return (abs(b2.ortho_beg - b1.ortho_end) == 1);
}




*/
